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ABSTRACT 


The periodic position and velocity perturbations of an 
artificial earth satellite are developed to the first order for all 
J , based on the theory by Brouwer as extended by Giacaglia. An 
explicit formulation is also provided for the subset J^, . 

The use of a position and velocity formulation circumvents the 
equatorial and circular orbit singularities found in conventional 
developments. The definition of the mean elements of the theory 
is modified to reduce the complexity of the position perturbations, 
as suggested by Merson's Theory, and the resulting changes to the 
secular terms are developed. In order to facilitate an empirical 
correction for drag, the observed mean motion is introduced as a 
mean element in place of the semi-major axis. 
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SECTION I 


INTRODUCTION 


The motion of a near-earth satellite is, in the first approximation, 
"Keplerian"; i.e., it conforms to certain empirical laws formulated in 
the seventeenth century by Kepler. In his Principia , Newton demonstrated 
that these laws described motion in an inverse-square force field. The 
force (or negative potential) function for such a field is of the form 

U » ^ (1) 

where r is the geocentric distance of the satellite and jj, is the 
gravitational constant. The motion is conventionally described by six 
"orbital elements," a, e„ I, M, uu, Q. Kepler's second law states that 
the motion occurs along an ellipse with one focus at the primary. The 
inclination, I, and the argument of the ascending node, Q, serve to 
locate the plane containing this ellipse. The eccentricity, e, and the 
argument of perigee, uu, define the shape of the ellipse and its 
orientation within the orbital plane. The semi-major axis, a, provides 
the scale of the ellipse as well as the orbital period; from Kepler's 
third law the period P is given by 





( 2 ) 


The location of the satellite within the ellipse is given by the mean 
anomaly, M, which measures the area swept out by the radius vector 
since perigee passage. In accordance with Kepler's first law the area 


1 




swept out and hence the mean anomaly increases at a uniform rate; then 


M(t) - M + n t 
o o 


where n is the mean motion, given by 


_ 2tt % -3/2 

n - — = y, a 
o P 


(3) 


(4) 


The mean anomaly must be converted to a geometric ang'.e to be of use; 
the eccentric and true anomalies, E and v, are related to M by 
Kepler's equation 


E - e sin E = M 

which must be solved by iteration, and by 


V 

tan - « 


l+e 

1-e 


, E 
tan — 


Following these computations, the geocentric position and velocity 
vectors _r and _r are given by (See Figure 1) , 

r = r U 
r = r U + rvV 


(5) 


( 6 ) 


(7) 

( 8 ) 


where 


r = 


r = 


l+e cos v 


% 


(—) e sin v 
P 


(9) 

( 10 ) 


• P* 

rv = (—) (l+e cos v) 
P 


(ID 
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POSITION IK THE ORBITAL PLANE 



POSITION AND VELOCITY IN GEOCENTRIC INERTIAL COORDINATES 


FIGURE 1 
3 









u 


cos u cos Q - sin u sin Q cos I 


( 12 ) 


cos u sin Q + sin u cos fi cos I 


sin u sin I 


V 


sin u cos 0 - cos u sin Q cos I 


(13) 


sin u sin Cl + cos u cos 0 cos I 


cos u sin I 


with 


2 


P m a?-| 


(14) 



(15) 


u = v-ht' 


(16) 


This simplified model cannot adequately represent the motion of a 
satellite; it must be augmented to include the effects of various 
"perturbing forces " These forces may be conservative, i.e , gravita¬ 
tional, or may affect the energy of the satellite's orbit, e.g., atmo¬ 
spheric drag and solar radiation pressure. More sophisticated models 
may be implemented via numeric integration or "special perturbations" 
techniques, which fall into three major categories: 


(a) Integration in cartesian coordinates of accelerations 
resulting from all forces acting on satellite, to obtain 
position and velocity (Cowell's method). 

(b) Integration in cartesian coordinates of accelerations 
resulting from perturbing forces only, to obtain deviations 
in position and velocity from a KepLerian orbit (Encke's 
method). 
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(c) Integration of the variations due to perturbing forces in 
the orbital elements of an "osculating 11 Keplerian orbit, 
i.e., the Keplerian orbit defined by the position and 
velocity of a satellite at each instant (Variation of 
Parameters). 

These methods have a theoretical disadvantage, in that the 
accumulation of roundoff and truncation errors must eventually result 
in inadequate precision; it appears that in practice the length of arc 
is limited by the uncertainty in modeling non-conservative forces, which 
involve not only a complex and highly variable atmospheric structure, 
but also the configuration and orientation of the satellite. A more 
practical problem in some applications is that the integration must be 
carried from epoch to the most distant observation, regardless of 
whether useful data exists in the intervening period. In addition, the 
integrated orbit, whether in terms of coordinates or elements, provides 
little insight into the effects of the various forces operating on the 
satellite, so that it is difficult to identify and correct deficiencies 
in the model. Despite these disadvantages, special perturbations 
programs are widely employed for precision tracking where the frequency 
of data mitigates their relative inefficiency, or the cost is justified 
by the requirements for maximum precision. They are also extensively 
employed in feasibility studies and similar investigations, to avoid 
time consuming (and possibly impractical) analytic developments. 
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It: is generally possible to obtain analytic expressions for the 
effects of the perturbing forces, to any desired precision and for any 
lime span Such "general perturbations" models are universally employed 
in routine cataloging systems, where a considerable number of satellites 
must be tracked with data that is sparsely distributed in time. In 
addition, general perturbations are usually employed in satellite 
geodesy. A considerable number of analytic theories have been developed 
for the conservative perturbing forces, i.e., the departure of the 
earth’s gravity field from an inverse-square law and lunar and solar 
gravity. 

For greater efficiency, semi-ana Lytic theories are often employed 
'or the luni-solar perturbations which are relatively small and of low 
requency; the same approach is generally followed for the solar 
adiation pressure perturbations. Thus, for solar radiation pressure, 
ehe analytic development may be carried through a formal integration of 
he perturbations, but the results are left as a function of the limits 
>f integration. These limits depend on the points at which the satellite 
enters and leaves the earth's shadow, which vary slowly with time. The 
evaluation of the perturbations proceeds n revolutions at a time, 

^ith the shadow limits reevaluated at each step. 

Analytic models of the drag perturbation have been produced for 
simplified atmospheric models The theory is complex, particularly 
when interactions with the earth’s oblateness perturbations are considered, 
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As a result, empirical models are generally employed, with only the long 
term effects of drag considered. The results are generally satisfactory 
for high altitude objects, but there appears to be considerable merit 
in the development of a semi-analytic drag theory. 


This paper deals only with perturbations due to the earth’s gravity. 
In geodesy the gravity field is described in terms of a reference 
ellipsoid, a reasonably tractable figure which approximates the figure 
of a rotating fluid in equilibrium to about 1 part in 10^. The actual 
gravity at any point is shown in terms of a map of the elevation or 
depression of the "geoid" with respect to this ellipsoid; this ”geoid n 
is an equipotential surface, i.e., a surface everywhere perpendicular 
to the local vertical. Before artificial satellites were launched the 
ellipsoid and "geoid 11 were determined from the reduction of direct 
gravity measurements and from astronomical determinations of the 
deviation of the local vertical from the local perpendicular to the 
ellipsoid This "geoid" data is not employed in tne theory of an 
artificial satellite, however. An analytic expansion for the potential 
is required; in spherical polar coordinates the generalized force (or 
negative potential) function is a series of Legendre polynomials and 
associated functions: 


0.4. s N 

r ' n=2 r 


n / 


n 

J P (sin 0) + E J P (sin 0) cos m (\-\ ) 

n n . nm nm nm 

m=l 


1 + 2 


n 

E 


n=2 m=0 


P (sin 3) (C cos + S sin m\) 
nm nm nm 


(17) 
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where is the equatorial radius and 3, > are the geocentric 

latitude and longitude of the satellite. The Legendre polynomials and 
associated functions are defined by 


P 

nm 


(x) 


2v m/ 2 ,n+m 

Slz*2 

2%! d x n4in 


(l-x 2 ) m/2 

2" j-0 


(- 1 ) 2 (2n-2 j ) ! x "' 1 ' 1 ' 2 
j! (n-j)! (n-m-2j)! 


(18) 


where I(~~) is the integer part of ( n, ~ ). The P (or P ) harmonics 
2 Z n no 

are "zonal," while the P harmonics are "tesseral " The largest 

nm 

-3 

coefficient is it is of order 10 . The remaining coefficients do 

not exceed the second order, i.e., 10 In order to evaluate the effects 

of the harmonics, it is necessary to substitute orbital parameters for 
r, {3, and In general, a method of successive approximations must 

be employed, so that a series of perturbations of increasingly higher 
order arises, e.g., 

first order: J 

2 

second order" J n , J , J 

2 n nm 

3 

third order* J n , J n J , J 0 J 

2 2 n 2 nm 

Most general perturbation theories neglect periodic effects of the 
second order, the residual perturbations will then be on the order of 
15 meters. However, under certain circumstances the perturbations due 
to higher order terms are amplified and must be included in a first 
order theory. If the potential function and its derivatives are 
expressed in terms of conventional orbital parameters, they will be 
found to have arguments of the form: 
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cos 


sin 


(n-2p+q)M + (n-2p)uu + m(Q-X) 


where n,m are the indices of the harmonic, p ranges from 0 to n 
(it is the parameter of a power series in sin I and cos I), and q 
ranges from -oo to °° (it is the parameter of a power series in e, 
with the lowest power of e being e^ ^ ) . From a simple first order 
theory, it will be found that M, uu, and Q all increase linearly 
with time, so that when the perturbations are integrated divisors will 
arise of the form 


(n-2p+q) n + (n-2p)(i + m(Q-X) 

o 

where n is the perturbed mean motion. The perturbations are classi- 
o 

fied in terms of n, m, p, q as follows: 


(a) Secular terms p = (n+q)/2 

q = 0 
m = 0 


These terms give rise to a linear increase in the elements 
M, uu, fi, and are therefore computed to the second order 
in a first order theory (so that the theory is valid for 
about 10-20 days, after which the neglected third order 
terms exceed the second order). These terms only arise for 
even order zonal harmonics, i.e., n=2, 4, ... ; the 

values of the even zonal harmonics are generally based upon 
observed secular perturbations. 
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(b) Long period terms p = (n+q)/2 

q t 0 

m - 0 

These terms have a divisor of the form -quo which is of 

-3 

order 10 ; second order forces therefore integrate into 

first order perturbations and must be included in a first 

order theory. There is no term of this form; if there 

were,a different type of solution would be required (there 

2 

is a term of this form which reduces to on 

integration). There is a special case for the “critical 
inclination" I & 63.4°, where 0) is of order 10 b , so 
that a "resonance" occurs In this case, either a special 
solution is employed or the long period terms are not 
integrated, i.e., they are left in the form of secular rates. 
The long period perturbations for even zonal harmonics are 
factored by the eccentricity e and can often be ignored; 
this is not the case for the odd zonals whose values are 
usually determined by analysis of observed long period 
variations in eccentricity and inclination. 

(c) Short period terms p (n+q)/2 

These terms have a divisor containing n^ so that the order 
of the perturbation remains unchanged upon integration. 
Therefore, only the terms need be included in a first 

order theory. 
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(d) Tesscral harmonic terms m 0 


There are two cases of interest here. For p = (n4-q)/2 
there are terms with frequencies near some multiple of the 
siderial rate, since to the zeroth order 
(n-2p)uj + m(Q-X) & -mX 

The integration results in an increase on the order of 

n /mX or about 16/m for near earth satellites. These 
o 

terms contribute perturbations on the order of 100 meters, 
and decrease in importance as n and ni increase. 

For 

(n-2p+q)n^ 4- (n-2p)uu 4* m(Q-X) « 0 

there is a resonance analagous to that holding near the 

critical inclination. The resonance will be in general 

larger for smaller values of (n-2p) and q. The principle 

resonances thus arise for 

m a* n n odd 

o 

where n is expressed in revolutions/day. Obvious cases 
o 

of potential near resonance are 24 and 12 hour satellites. 
High order resonances, e.g., (n,m) of (13,13), (15,13), 

and (15,14) have been reported for certain satellites with 
magnitudes on the order of 100-150 meters and periods of 
2.5-5 days Obviously, by going to a sufficiently high 
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order harmonic, a resonance can be found for any near 
earth satellite. Fortunately, the net effect of these 
higher order terms is considered to approach the second 
order. 

This paper does not deal with the tesseral harmonic perturbations. 

It is limited to some minor modifications of the secular terms developed 

by Brouwer^) and extended by Giacaglia^^ , and to a non-singular 

development of the long and short period perturbations due to the zonal 

harmonics. In the Brouwer and Giacaglia papers the perturbations of 

conventional elements are computed, which leads to singularities for 

(3) 

low eccentricity or inclination Lyddane showed that the problem 

could be circumvented by either computing the perturbations to "non- 

singular" elements, e.g., e cos M and e sin M, or by computing the 

perturbations in the position and velocity vectors. The former approach 

(4) 

is employed in most general perturbations ephemeris generators , while 
the latter approach is employed in this paper. Although Garfinkel^^, 
Kozai^\ and Merson^^ have computed some of the position perturbations, 
velocity perturbations have generally been neglected. 

The use of position and velocity perturbations has the advantage 
of revealing the "real" or observable effects of the perturbing forces; 
Merson^^, for example, has shown that some of the apparent perturbations 
of the orbital plane affect only the velocity vector and can be ignored 
in a tracking network based on positional data. In addition, a 
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position and velocity theory appears to be somewhat more efficient than 
a "non-singular" elements theory, particularly when only positional 
data is used for element correction. The position and velocity theory 
has one disadvantage, in that the frequency of the long period terms 
becomes comparable to the short period terms, and they must be recom¬ 
puted for each ephemeris point. (However, they are recomputed for 
each point in most theories, whether or not the computation is necessary.) 
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SECTION II 


PERTURBATIONS IN POSITION AND VELOCITY 


In developing the perturbations it is convenient to use the 
momentum unit vector W, given by 


W = U x V = 


sin 0 sin I 
- cos 0 sin I 
cos I 


The perturbed position and velocity may be computed as 


r = (r + 6r) (U + 6U) 

r = (r 4- 6r) (U + 6U) 4- (rv 4- 6rv) (V 4- 6V) 
or the perturbations alone may be computed as 
6r = 6r U 4- r 6U 

6r = 6rU + 6rvV+r 6U + rv 6V 
ignoring second order terms. 


The quantities 6U and 6V may be written as 

6U = V (6u + cog I 60) + W (sin u 61 - cos u sin I 60) 

6V = “JJ (6u + cos I 60) + W (cos u 61 + sin u sin I 60) 

and hence we have 

6r_ = 6r U + r (6u + cos I 60) V 

4* r (sin u 61 - cos u sin I 60) W 

6r_ = | 6r - rv (6u 4- cos I 60) J U 

H-16 rv + r (6u 4- cos I 60)J V 

4- j r (sin u 61 - cos u sin I 60) 

4“ rv (cos u 61 4- sin u sin I 60) W 


angular 


(19) 

( 20 ) 
( 21 ) 

( 22 ) 

(23) 

(24) 

(25) 

(26) 


(27) 
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For 6r, &r, and 6u we can either use Taylor series expansions 

/g\ 

in the conventional elements, or the ingenious equations of Izsak J with 
the Brouwer determining function S written in Hill's canonical 
variables {r, G, H I r, u, 0} . The equations are 


c Sa c , a sin v 

or = r — - a cos v oe + - e oM 

a T| 


(28) 


6 r = 


• a ,, 1/2 v 2 ,, 1/2 .2 

r 6a . N (1+e cos v) . , . v (1+e cos v) 

- -x — + (-) -^- *— sin v 6e + (-) cos v e 

2 a P 2 P 3 

iq ri 


6u = 


(2+e cos v) c (1+e cos v) . s ,., 

^sm v Oe + --- OM + OUD 

T) *1 


(29) 

(30) 


or 


6r = 
6 r = 
6u = 


9S 

5r 

3S 

5r 

3S 

5G 


(31) 

(32) 

(33) 


The second set of equations appears much simpler, and has been 
solved for the short period terms by Izsak. They are not so simple, 
however, when dealing with the long-period terms containing trigonometric 
functions of 'JJ. 


The perturbation in rv may be computed from 


o rv = rv 


_1 5_a cos v (1+e cos v) - 


2 a 


6 e 


T1 


sin v (1+e cos v) < w 
- 2- e 6M 

T1 


(34) 


Brouwer and Giacaglia employ the Delauney canonical variables 
(L, G, H j 1, g, h], where 

L = (p.a) 2 1 = M 
G = LT) g = 0J 
H = G cos I h = 0 


5M 
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or 

£ £,G, . / 6G 6 Tv /OC\ 

6rv = 6(-) = rv (—--) (35) 

r u r 

using intermediate results for 6G. 


If the Taylor series expansions are employed, it is possible to 
rewrite Equation (27) using Equations (29), (30), and (34) as 



1/2 


x 


e sin v 
2 


6 a 
a 


(1+e cos v) sin v P 
2 C 

n 



' ^~^ ~ 3° S " (1+e cos v) (6uu + cos I 60) 

r\ 


1+e cos v 6a . cos v+e c . . /P . . T e ns 

-_-j. --- 6e + e sm v (6uu + cos I 6U) 

z a z 

T| 

(cos u + e cos uo) 61 + (sin u + e sin uo) sin I 60 


U 


V 



(36) 
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SECTION III 


SHORT PERIOD PERTURBATIONS 


/O') 

Izsak v has already computed 6r, 6r, and 6u as 


j rs a l 

6r = ssin I cos 2u + (1-30 ) 


1 + 


M 


e cos v 


1+e cos v 1+n 


(37) 


6r = - [i 


* J 2 a e 


4 p 


5/2 


6 u = 


J 2 a e 

8 p 2 


2 2 
2 sin I (1+e cos v) sin 2u 


. /, ^^2 n . , (1+e cos v) 2 \ 

+ (1-39 ) e sin v I 0 + ^ > 


9 9 1 «. qq 

6(1-50 ) (v-M) + 4(1-60 + •) e sin v 


(38) 


+ (1-30 ) (1-0) sin 2v + 2 (50 - 2) e sin (2u-v) 


+ (70 2 - 1) sin 2u + 20 2 e sin (2u+v) ) (39) 


where 


- coefficient of the second zonal harmonic 


earth's equatorial radius 
9 = cos I 


a = 
e 


From Brouwer's theory, with 

J 2 a 2 
y = 2 e 

2 2 a 2 
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we have 


5 r ^? a e 2 / e 

— = 3 - — sin I ( cos 2 u + e cos ( 2 u-v) + — cos ( 2 u+v) 

G 4p ' 3 


(40) 


^ O a p ° 

cos I 6 Q = - -— 0^ \6 (v-M 4- e sin v) - 3 sin 2u 

4 p Z 


- 3 e sin (2u-v) - e sin (2u+v) 


(41) 


61 = ■ — e - sin 19(3 cos 2u + 3 e cos (2u-v) + e cos (2u-K;) 

4 p 


Hence 


(42) 


6 u + cos I 60 = 


^ 6 \sin^ I) | + 2 e sin ( 2 u-v) 

4 p‘ ■ 


(1-39 2 ) 


2+rv 

3 (v-M) + 2 e sin v (y—) 


+ -^2 sin 2 v 


(43) 


J 2 a e 

sin u 61 - cos u sin I 6 Cl = -— sin I 9 

4 p 2 


X 


- 3 sin u - 4 e sin <ju 


+ 4 e cos u sin v + 6 cos u (v-M) 
J„a 2 


(44) 


cos u 61 + sin u sin I 60 = ———-j - sin 19(3 cos u + 4e cos (u+v) 

4 p 2 


- 6 sin u (v-M) ) (45) 
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6 rv = M* 


h J 2 a e 


4 p 


5/2 


(1+e cos v) 


X 


(sin^I) 


- (1-30 2 ) 


2 cos 2u + 2e cos (2u-v) + e cos 2u cos 


•| (1+Tl) + e cos v cos 2v 


(46) 


Substituting these terms in Equations (26) and (27) gives, after 
some simplification, the short period terms in Table I. 

In some tracking programs the relationship of the mean semi-major 
axis and the secularly perturbed mean motion is taken from Kozai's 
equation 14^^: 


J„a 


^_3 .. f, . „ 2 e „ _2 


7T = \i |1 + 3 


(1-30 ) T) 


4 p 


Since 


n = n 


this implies 


(47) 


^2 3 e 2 

1-3 - ■ - (1-30 ) T| + . 

2 P 


(48) 


a = a 


1 + 3 


J o a 
2 e 

,p 2 


(1-30 ) T1 + . . . 


(49) 


where a is the mean semi-major axis of the Brouwer theory, defined by 

n 2 a 3 = M> (50) 

o 
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SHORT PERIOD PERTURBATIONS 
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The use of a in place of a in the computation of unperturbed _r 
and jo requires that a compensating perturbation be applied to 6jr 
and 6jr: 

^ - - 3 (1-39 2 ) n (51) 

a a , l 

4 p 

This results in the following changes in Table I: 

(a) The term 2T) in the U component of 6jr 
becomes - T\ 

(b) The term -T] e sin v in the U component of $£ 

T\ 

becomes 4- — e sin v 

Merson has developed a theory^^ in which the short period position 
perturbations are minimized. His formulae (153-157) relate the osculating 
elements to conditions at the ascending node, to the second order in J^. 
By eliminating all first order terms whose argument is a multiple of uu, 
and terms factored by u, a set of first order pseudo short period terms 
is obtained. These result in the position perturbations given in 
Table II. (The first order terms factored by u in Merson's theory are 
actually the sum of secular terms factored by M and short period terms 
factored by (u-M). They can therefore be included in the secular term 
computation, and we have taken this approach.) 

Now, if we define £ . to be Brouwer’s mean elements updated for 

i 

secular perturbations, and to be "smoothed" mean elements updated 
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SHORT PERIOD PERTURBATIONS FOR "SMOOTHED" ELEMENTS 


>1 31 
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for secular terms, then: 


r(e j) + 6r 


Brouwer 


i) + ^ (e + 6r, 

i-1 Gi 


i 1 —Smoothed (52) 


so that we have three simultaneous equations: 


V ?- Ac. - 


i=l 


ae i 1 


6_r - 6 _r 

Brouwer Smoothed 


(53) 


where 


Ae . = e ! - e . 

i i i 


Using Equations (26), (28), (30) and Tables I and II, 


Aa * . a sin v 4W | „ ,, ^2 X V 

r — - a cos v Ae H-eAM = r (X- (1-39 ) A 

a rill 


2 

l+3ri+e 

+ (2+n' 

1 e 2 

e cos v + — cos 2v 


CM 

_1 

i 1+n J 

1 2 



. e . „ \ Ae . (1+e cos v) 
2 sin v + — sin 2vJ — H-^- 

*n T1 


Am + Ago + 9Afi 


= - r 


a 1 (l-30 


r- 

i 

/ 2 

“ 

3 (v-M) + 

2+n 

1 e 

2e sin v + — sin 2v 



i 1+ n i 

l 2 1 

- 


r (sin u Al - cos u sin IAQ ) = r Ot sin I© 


6 cos u (v-M) - 3 sin u 


( 54) 
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If the Ae are restricted to contain only constants or terms with 
(v-M) as angular arguments these equations may be solved to yield: 


- = 2a. (1-39 2 ) T) 
a l 


A e = 

- a. 

L U-39 2 ) 

2+t] 

i+n 

el] 

AI = 

- 3 

a sin 19 



AM = 

0 




Au) = 

- 3 

a i (1-59 2 ) (v- 

■M) 

AQ = 

- 6 

a 9 (v-M) 



The resulting changes in the velocity perturbations may be obtained 
from an equation similar to Equation (53); the signs of the Ae . 
should be changed to yield: 


6jr - 6r 

Smoothed Brouwer 


i 

i=l 


Sr . 
r* Ac, 
9e i i 


(56) 


With the aid of Equation (36) 


6r - 6 £ 

Smoothed Brouwer 


<h X 


(1-30") 


e sin v 


■T1 + 


2+Q 

i+n 


(1+e cos v) 


+ 3(l+e cos v)(v-M) 


+ (1-39 ) 


- 1 - e 


e cos v 


T1 + -^j|- 3e sin v (v-M) 


- 3 sin 19 


cos u+e cos UL) + 2(sin u+e sin uu) (v-M) 


U 

V 

;■ 


( 57 ) 
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which, when added to 6_r from Table I, yields the velocity perturbations 
in Table II. 

These position and velocity perturbations represent a useful 

simplification of the equations in Table I; corresponding modifications 

2 

to the and secular terms are derived in Section V. 
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SECTION IV 


LONG PERIOD PERTURBATIONS 


The first order long period perturbations for J^, n> 2, have 

(2) (9) 

been developed by Giacaglia and Garfinkel and McAllister 

Giacaglia's results have been employed since they are more readily 

truncated; the published paper unfortunately contains a number of errors: 

(a) For the even harmonics, the upper index limit for k in 

Z is incorrect, since values of i within the permitted 
kij p 

range do not exist for k= Furthermore, for k = 0, i 

can exceed j; in this case, 6 ,^vanishes.* 

p+l, 2j, 2i 

Therefore 


£ 

kij 


2z2 

2 

£ 

k=0 


min 

£ 

i=l 


-2k ’ 


Ell 

2 

£ 

j=l 


Similarly for the odd harmonics 



szl 

2 

min 

UJ 


£ 

- £ 

£ 

l 2 1 

f £ 

kij 

k=0 

i=0 


j=0 

where 

P-2.1-1 

2 

was 

incorrectly ; 

1 imit 

of i 





These practical problems reflect theoretical errors; i.e., k = -^ terms 
are secular terms aid also appear in the £ portion of the disturbing 
function. kj 
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(b) In equation 21, the term 
„ 1-e 2 


should be 


- 2j 


, g-o . 


(c) Equation 22 must be multiplied by H. 

(d) In equation 28, the factor 

lizi 

2i+l 

should be 

2j+l 

2i+l 

(e) In equation 30 the term 

p-2k 

• 2 T 

sin I 

should have a minus sign, and the factor 

cos (2i+l)g 

21+1 

has been omitted. 


Also, Giacaglia has the wrong signs for the third through fifth 

harmonic Brouwer coefficients (k , v , and A -) as functions of J . 

n T n n. 0 n 

With these corrections, Giacaglia's results may be summarized as 

6e - s Ae (58) 

where S and Ae for even and odd n are summarized in Table III. 
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OIACAC LIA'S LONG PERIOD PERTURBATIONS 
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co 
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It is convenient to reduce these equations to a common form. First, 


note that the summations < 

can be written 

as 


n-2 n-2 

n-2i 



2 2 

2 



£ £ 

Z 

n 

even 

i=l j=i 

k=0 



n-3 n-3 

n-2i-l 



|cs 

2 



z z 

£ 

n 

odd 

i=0 j=i 

k=0 



Now introduce a variable \ 



X = 2i 


n 

even 

= 2i+l 


n 

odd 

so that 




II 

ro 

• , (n-2) 

n 

even 

= 1, 3, 

, (n-2) 

n 

odd 

Similarly introduce 




u = 2j 


n 

even 

= 2 j+1 


n 

odd 

so that 




p. “ \, ^+2 

, ... , (n-2) 



Also introduce 




V = n-2k 





(59) 


(60) 


v = X , X +2, . 


n 


( 61 ) 




Note that S is factored by e sin I. Moving this factor from S 
into the As, and employing 

% = n/2 - id (62) 

we obtain the results in Table IV. 


It is now a simple matter to obtain 


6 r = - aT] sin IS 


{ 


COS V cos 


+ \ 


sin v sin 


«} 


(63) 


c . (JL 1/2 2 f p. l 

°r = (—) (1+e cos v) sin I S «j sin v cos X§ — cos v sin X§ | ( 


64) 


6rv = (~)^^(l+e cos v) sin I S 


(1+e cos v) 


cos v cos X§ + r- sin v sin X§ 

X 


- e cos X 
(65) 




6u + cos I 6Q = sin I S < (2 sin v+e ) cos X? 


. cos 2v . 3e N P* sin X? 

- (2 cos v+e — - + — ) ---“ L 


sin u 6 I - cos u sin I 6 fi - -e 9 S 


+ e (2n-5) } 


( 66 ) 


X sin u cos X? + 


cos u 


v - 


10 sin I 
2 

1-50 


in X? 


sin 


} 


cos 


u 6 I + sin u sin I 6 ft = -e 9 S 


(67) 


*{ cos u cos X? + sin u 


10 sin I 


- v 


sin Xg | 


1-50 

from which, with Equations (26) and (27), the results in 
Table V are obtained. 


( 68 ) 
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MODIFIED LONG PERIOD PERTURBATIONS 
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LONG PERIOD PERTURBATIONS IN THE COORDINATES 


>1 


>1 


^1 


^1 


^1 




H 

C 

•H 

CO 

+ 


•M 
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It will be noted that S vanishes for n=2. 'll* this were not the 

case, it vjould be necessary to modify the basic unperturbed solution, 

since the resulting terms would have a J ^/J coefficient and could not 

be treated as a perturbation. However, the magnitude of J is such 

that the long period solution for must be carried to one higher 

approximation than for the other J^, yielding perturbations with the 
2 

coefficient J^/ J 0 or . As giver: by Brouwer, these terms are: 


6 a 

6 e 

61 

fiM 


* 0 

2 2 

= r\ sin I cos 2uu 
2 

= -S 2 e 9 cos 2ju 
o 3 

. o 2 11 . T . , 

- — sin I sin 2o) 

z e 


(69) 

(70) 

(71) 

(72) 


6uu 


S 2 6 


sin I 




s in I 


_ 109 2 1 

(1-59 2 ) (1-159 2 )J 


in 2uu 


(73) 




e sin 19 


|_J_ + -12- 

1 sin 2 I (1-5S 2 ) (1-153 2 ) 


sin 2ju 


(74) 


where 


S 


2 

2 


1_ 

16 


2 

P 


e s in I 
(1-59 2 ) 


(1-159 2 ) 


(75) 


The additional perturbations are easily obtained: 
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6 r„ = r -77 ————- ———(1-150 2 ) ( - sin I (1-t-e cos v) cos (2u-v) U 
" 2 10 p 2 (1-59 Z ) 


-!- sin I 


2 sin (2u-v) 4- 


c sin 2 u 


e9 


sin (u-2id) - 


10 sin""I cos u sin 2i) 

(1-59 2 ) (1-159 2 ) 


W 


(76) 


^—2 = 


d-159 2 ) 

?! 16 p 2 (1-50 2 ) 


X 


- sin I (1+e cosv) 


sin ( 2 u-v) - 


e sin 2 a) 


-I- s in I 


cos ( 2 u-v) + e cos 2 U) - 


2 

e sin v sin 2'JD 


- e9 


cos (u-2uu) + e cos 4- 


10 sin I (sin u+e sin w) sin 2 uu 

(1-59 2 ) (1-150 2 ) 


W 


(77) 


In most general perturbations formulations it is customary to employ 
2 

only the J^ 9 anc * ^4 lor.g period perturbations. For comparative 

purposes, these are given explicitly in Table VI. An evaluation yields 

For J 3 


X * P* = 1 V = 1,3 


, a -1 

S 3 - + 2 37T (1 - 59 > X 


-5 sin I 



(78) 
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a 

si 3 e 

6 r 0 = r- -- 

“3 2J p 


sin I (1+e cos v) sin u U 


+ sin I (2+e cos v) cos u V 
+ e9 cos v W 


( 79 ) 


»i 3 - 



\ J 0 a I 

3 e ) 

\ p 1 

ai 2 P \ 


sir I (1+e cos v) cos u U 

- sin I (sin u+e sin uu) V 
- e9 sin v W 


(80) 


For J, 


X = H - 2 v = 2,4 


S = + — e - sin 1 V l 2 

4 16 J 9 p 2 (1-50 2 ) ^ 


6 v = 2 

7 sin*" I V = 4 


(81) 


- J, a . _ 

t 5 4 e e sm I 

or = r — - - 

- 4 16 jV d-50 2 ) 


(1-79 ) < - sin I (1+e cos v) cos (2u-v) U 


+ s in I 


2 sin (2u-v) + 


e sin 2u 


V 


- e9 i sin (u-2^) 


2 sin I cos u sin 2ju 


(1-50 2 ) (1-70 2 ) / 


W 


(82) 
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^—4 = 


if* JL J 4 a e e sin I ,, ?fl 2. 
P l 16 J, P 2 d-59 2 ) U } 


^ < -sin I (1+e cos v) 


sin (2u-v) 


e sin 2uu 


+ sin I 


cos (2u-v) + e cos 2UJ - 


sin v sin 2uu 


eQ cos (u-2^ 1 ) + e cos UU + 


2 sin I (sin u+e sin W) sin 2W 

(1-59 2 ) (1-79 2 ) 


W 


(83) 


From which the equations in Table VI are easily obtained. 
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LONG PERIOD PERTURBATIONS 
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SECTION V 


SECULAR TERMS 


The results of Brouwer^ and Giacaglia^^ are given in Table VII, 


where 


n-2 

. X* r~ (-l) k (2n-2k)! e 2ji sin"' 2k I 




n-2 

n even 


P “ 0 2 2(j_k) k:(n-k):(^!) 2 (j!) 2 (n-2j-l)! 


(84) 


In most practical applications it is the parameter n rather than 
a which is determined. From the relationship 


1 i . 2 , 

ua = 1 ' a +a + ••• 

where (Ha) represents the terms factored by n Q in the equation for 
"n we obtain to the second order: 


n = n (1-a + a +...) 

o 

- * < 1 + ’ 11 (1 - 3 « 2) * tM ^ 71 ' 2 


4 P 2 

, J 2 . 4 

3 2 I-Tl 


3-309 2 + 359 4 


128 


10-8T)-25e 2 + (-60+48tl+90e 2 )e 2 + (130-72n-25e 2 )e‘ 


1 1 - St| (2J , - 3) + 

j- term above j ) 

l 

L J J 


( 85 ) 
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SECULAR TERMS ACCORDING TO BROUWER AND GIACAGLIA 


i 



39 


term above 
































if n is substituted icr in the equations for 'JJ and Cl in 

2 

Table Vll. fit- eo. j L of the terms is reduced; they become: 


Tit 


T 2 * r 
T g 

n jig '~~1T '- 1() - 25eJ + (- 36 + 126 e 2 ) 0 2 + ( 430 - 45 e 2 ) 9 4 


(86) 


For O 


r .0 


.2 4 f 

_3 J 2 d c 

32 4 

P 


> 9 


4 - 9 e + (- 40 + 5 e^) 9 “j 


( 87 ) 


• b. competed from the relationship 


a - 2/3 =a 1/3 r,- 2/3 (1+2/3 a - 1/9 a 2 ...) 

o 


1)1 * ..h vc 


a - 


1/3 

U r. 


-2/^ j ^2 a e 2 

< 1 - 1/2 ri ( 1-39 ) 


+ — 


1 J 2 a 6 


^4 p4 


r\ 


10 + 12 n- 25 e 2 + (- 60 - 72 Tl+ 90 e 2 ) 9 2 
+ ( 13 O+lO 8 ri- 25 e 2 ) 0 ' + 


15 J 4 d e „ 2 

64 — 71 6 


3 - 3 O 0 2 + 359 ^ 


V , 173 n 


- 2/3 


i+s"n 


2 \ 

. 2 

3 e 2 


term above 


( 88 ) 
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Note that p, which is a function of a, enters into these equations. 
If we employ 


1/3 —-2/3 2. 

p = p, n (1-e ) 


= P 




T1 (1-39 2 ) + .. 



in the equations in Table VII five J T\ terms are added to Equations 

(86), (87), and (88). It is therefore more efficient to compute p 

-2 

from n and p from 


-2 

P 


-2 


J 2 a e 2 

1 + T| (1-30*) + ... 


(89) 


~ 2 

for use in the terms. Either p or p may be used in the J 0 

terms, since the error in using p is of the third order. If, however, 


Equation (88) is only carried to the first order, it may be entered with 


p, and p computed from a, so that Equation (89) is not required. 


If the Kozai a is to be employed as an element, then Equation (88) 
is replaced by 


a - a (1 + T 


3 J 2 a a 


4 2 

P 


T1 (1-39 2 )) 


1/3 -' 2/3 

|i n 


J 2 a e 2 

1 + T) (l - 39 ; 

4 p 


(38a) 


41 








t » the first order. The use of n and 


P 


a (1-e ) 


in the equations in Table VII would again lead to additional terras 
in T| in Equations (86) and (87). The preferred approach is therefore 
to correct ~p using 


-9 

P 




2 ^ 

(1-39 ) 


(89a) 


A rather more complicated approach is employed in the GEPERS 
(4) 

luatiue developed for use in the SPACETRACK system. The parameters 


tiylcyed include a and p, with a mean motion parameter n defined 

t - 

2 


na -u 1 T| (1-30U 

P 


(90) 


t./ich is interpreted to mean 


n = n 


, 3 J 2 a e 2. 

1 • 4 —~ in a- 30 ) 

p 


= n 


1 ’ 4 11 (1_302) 


T 2 4 

o o / 

+ Y^g T) (48 T 1-288T|0 Z + 432T10 4 ) 

P 
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(where either p or p may be used in terms of order J^.) ^ 

this equation is subtracted from the equation for “n in Table VII 

and n is subtracted for n in the second order terms we obtain: 

o 


- - n ' 1 + 


10-3Zn-25e 2 + (-60+192Tl+90e 2 )9 2 
+ (130-288r|-25e 2 )9 4 


J a 4 

45 4 e _ 2 

128 _4 r ' 6 


3-309 2 + 359 4 


(91) 


The expressions for uij and Q in GEPERS employ n and "p in 
place of rT and p; "n agrees with H to the first order, but 
the use of "p makes it necessary to modify Equations ( 86 ) and (87). 
The results may be obtained from Equation ( 86 a) and are: 

For jj 


t 2 * 

3 J 2 a e 
128 p 4 


• 10-48H-25e 2 + (-36+384n+126e 2 )9 2 + (430-720n-45e 2 )9 * 


( 86 a) 


For Q 


n9 


,2 4 

3 J 2 a e 


2 . „2 


| 4-24n-9e + (-40+72fl+5e )9 


(87a) 
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Suppose now that we wish to use the smoothed elements n, a', 
e', I', 'j ) 1 and Q' as defined in Equation (55). If we enter the J 

4 

terms for uj and Q in Table VII with 
P' - a' (1-e 12 ) 

V = (i-e 12 > 1/2 

0 ' = cos I' 


then the terms in Equation ( 86 ) and (87) will be changed by 

2 


6 uu = n 


3 ^2 a e 2 

{ (1-59 Z ) 


^ 6 a . 4e 6 e , 10 sin 10 61 
2 _ + __ + 


r\ 


(1 - 50 ) 


(92) 


60 = n 6 


3 J 2 a e 2 | 


0 6a . 4e 6e 

sin I 61 

2 2 
\ p i 


- ^ + o 

a i 2 

L 1-e 

e 


(93) 


where the £e^ are related to the Ae of Equation (55) by 


6 e . = e . - e ! = - A£ . 
111 1 


The results are 

For uo 


n 


3. 

128 


T 2 4 

J 2 a e 


4 


P 


-42-57e 2 + (-20+382e 2 )9 2 + (190-525e 2 )6 4 


( 86 b) 
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For Q 


t 2 4 

-■ 3 J 2 a e 

n9 32 4 


2 2 2 
■25e + (-4+53e )9 


(87b) 


These terms differ from those given by Merson^^ because he employs 
u as the argument of the secular terms rather than M (in the form 

nt ), i.e.j his rates apply for a nodal rather than an anomalistic 

2 

* ^ yjQ 

frequency and may be obtained by subtracting — and — from 
Equations (86b) and (87b). 


Now since from Equation (55), 

2 


J a 

UU ' - UJ = - 7 -Mr- (1-59 2 ) (v-M) w— (v-M) 
4 Z n 


J a * 

Q' - n = - | 9 (v-M) (v-M) 


if we define the increase in mean anomaly to be 


AM = n At 


(94) 


and we define Av to be the related quantity 


Av = AM + (v-M) (95) 

(note that Av is not zero at the epoch, unless epoch is defined to 
be at perigee or apogee ) 
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we have 


uu' = 


• uu / v 

uu + uu At H—=“ (v-M) 

o n 


(96) 


uu 1 +^=- Av 
o n 


and similarly 


Q ' = 0' +4 Av 
o n 


(97) 


where and — are obtained from Table VII as modified by 
n n J 

Equations (86b) and (87b), It remains to determine a* 

(and, thereby p 1 ). Since 

i J ? a 2 
1 2 e 


a' = a 1 + - — f-T) (1-39 ) 
1 P 


Equation (88) becomes 


, 1/3 -- 2 « J, . i J k‘ 

‘ " < l +64 -T 11 

P 


10-4Tl-25e 2 + (-60+24n+90e 2 )9 2 


+ (13O-36r|-25e 2 )0 4 


15 J ** e r. 2 

4 ^ 6 

P 


3-30Q 2 + 350 4 


1/3 -' 2 / 3 

p. n < 1 + ST| 


4j 1-e 

3 2 

J e 


- 2 + 


J term above 


(88b) 


(n > 2) 
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Thus, to the first order, a 1 may be computed from n using the 
Keplerian relationship and no special provisions are required for 
computing p*. 

The non-conservative perturbations are sometimes represented 
by polynomials in time. If, for example 


• ry •• o 

M = M + n At + r At + y At 
2 6 


(98) 



(99) 


then the coefficients for a 



( 100 ) 


are easily obtained from Equation (88b) as 



( 101 ) 


a 



( 102 ) 


d _ -J i.1 -L 

2 " i 9 I n " 3 n 3 


(103) 


A similar polynomial can be employed for e 



( 104 ) 
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The time dependent terms may be obtained empirically, or from the 
approximation that perigee height remains constant 


d- 6 ') -° 


which yields 


e = 


1-e 

a 


a = - — 


2 1-e 


(105) 


2 


1-e 

a 


a 

i 2 


a 

a 


1-e 

n 


1 * 2 
1 n 


1 .. 

3 n 


(106) 


The inclination is usually assumed constant, although it may be 

represented by a similar series. Where these secular perturbations 

on n, a 1 , and e* exist, they must be reflected in the motion 

of node and perigee. The variations in TT are already incorporated 

in Equations (96) and (97). It is usually adequate to carry the 

additional variations due to a 1 , e 1 , and I* to the first order 

in J * the necessary derivatives of —=r- and -O r- may be obtained 
2 n n 

from Equations (92) and (93). . The resulting equations are 




f \ 


a 


U) 





o 


n 


as 


+ 

4 

n' 


n 


a 


o 


h 

k , 



- 



/ 10 sinlQ^ 


Av 1 


„ 4e 

2 — + _— 

a i 2 

1-e 


£e + 


l-59 z 

sin I 

9 


si 


(107) 


where the 6's represent the change from epoch, or 
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w 


n 1 


J 


uu 


i; 


o 


+ 


4L 

n 


n 


Av (1+a) 


(108) 


where (using Equations (102) and (103) ) 


a = 


4 n . 4e 

- - + -r e + 

3 n .2 
1-e 




/ \ 

10 sip 19 

1-50* 

sin I 


0 


\ 


2 n . 2e ... 

+ -r - + -~ e + 

3 n . 2 

1-e 


i At 


^ 10 sin 10^ 
1-50 2 


sin I 

0 


\ 


.. \ 

1 

2 


At 


/ 


(109) 


or, for constant perigee height and inclination (see Equations 
(105) and (106) ), 


a 


R A t + ( - 1Q - 6e 

3(1+e) n C 9(l+e) 


R \ 2 + lilzRl r\ 

nI 3(1+e) n j 



(109a) 


In some cases it is desirable to express the secular terms as 
functions of time rather than as functions of Av. It is possible 
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to expand Equations (108), (109) and (109a) using Equation (95) 
and by separating like powers of At to obtain 



At + 


v-M 


/CA 

UJ 

2 


u 

2 

\ 


At 2 + 


UJ 

6 


Q 

K 6 * 


At 3 + . 


where 


/.\ 


A\ 

UJ. 1 




UJ 








T\ 





s 

• 

n 



n 


£L 




/ 


n 






> 






/ \ 


/ . \ 

UJ 


iL 


10 sin 10 


2 


n 
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